/*
 * !---------------------------------------------------------------------
 * ! Random number generator functions for Kruununhaka 
 * !---------------------------------------------------------------------
 * ! Copyright (C) 2010 Jussi Lehtola, University of Helsinki
 * !
 * ! This program is distributed under the terms of the GNU General
 * ! Public License as published by the Free Software Foundation;
 * ! either version 2 of the License, or (at your option) any later version.
 * !
 * ! This program is distributed in the hope that it will be useful,
 * ! but WITHOUT ANY WARRANTY; without even the implied warranty of
 * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. In addition, the 
 * ! author must be properly cited whenever the program or some parts 
 * ! originated on it are used. 
 * !---------------------------------------------------------------------
 */

#include <math.h>
#include <gsl/gsl_rng.h>

/* Random number generator */
gsl_rng *gen;


/* Allocate RNG */
void rnd_alloc_(void) {
  /* Use Mersenne twister generator */
  gen=gsl_rng_alloc(gsl_rng_mt19937);
}

/* Deallocate RNG */
void rnd_free_(void) {
  gsl_rng_free(gen);
}

/* Set RNG seed */
void rnd_set_seed_(int * seed) {
  gsl_rng_set(gen,*seed);
}

/* Get uniform number */
double rnd_uniform_(void) {
  return gsl_rng_uniform(gen);
}

/* Get Gaussian random number using Box-Muller algorithm */
double rnd_gaussian_(void) {
  /* Numbers to generate within square */
  double v1, v2;
  /* Distance from origin */
  double w;
  /* Radius */
  double r;

  /* Variables for second random number */
  static int saved=0;
  static double gnum=0;

  /* Have we got a number in store? */
  if(saved) {
    saved=0;
    return gnum;
  }

  /* Generator */
  do {
    /* Generate numbers within square [-1,1]x[-1,1] */
    v1=2.0*gsl_rng_uniform(gen)-1.0;
    v2=2.0*gsl_rng_uniform(gen)-1.0;
  
    /* Check if they are within the unit circle */
    w=v1*v1+v2*v2;
  } while(w>1);

  /* Calculate radius for Gaussian numbers */
  r = sqrt(-2*log(w));
  /* The wanted random numbers are */
  saved=1;
  gnum=r*v2/sqrt(w);

  return r*v1/sqrt(w);
}
